#include "current.h"
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf_gamma.h>
#include <assert.h>

/**
 * Calculates \Phi_\omega shape function Fourier spectrum 
 * @param shape - shape object to study
 * @param w - radial frequency
 * @return spectrum value at the specified frequency
 */
static complex shape_get_spectrum(const union current_shape_u * shape, double w){
	complex ans = GSL_NAN + GSL_NAN*I;
	switch(shape -> id){
		case CURRENT_SHAPE_EXPEXP:{
			double Q = shape -> expexp.Q;
			double alpha = shape -> expexp.alpha;
			double beta = shape -> expexp.beta;
			ans = (alpha*beta*Q) / ( 2*M_PI * (beta-I*w) * (alpha-I*w) );
		}
		break;
		case CURRENT_SHAPE_POWEXP:{
			double Q = shape -> powexp.Q;
			double alpha = shape -> powexp.alpha;
			uint16_t n = shape -> powexp.n;
			ans = Q / (2*M_PI) * cpow(alpha/(alpha-I*w), n+1);
		}
		break;
		default:
			assert(FALSE);
	}
	return ans;
}

/**
 * Calculates tau value, meaning signal start delay
 * @param delay - object to study
 * @param l - length to calculate delay for
 * @return tau
 */
static double delay_get_tau(const union current_delay_u * delay, double l){
	double ans = GSL_NAN;
	switch(delay -> id){
		case CURRENT_DELAY_LINEAR:{
			double v0 = delay -> linear.v0;
			ans = l / v0;
		}
		break;
		case CURRENT_DELAY_EXPONENTIAL:{
			double v0 = delay -> exponential.v0;
			double gamma = delay -> exponential.gamma;
			ans = -1/gamma*log(1 - l*gamma/v0);// If L is too long, it will produce NAN
		}
		break;
		default:
			assert(FALSE);
	}
	return ans;
}

/**
 * Calculates amplitude decrease coefficient ksi.
 * @param decrease - object to calculate decrease value by
 * @param l - length to calculate value for
 * @return ksi
 */
static double decrease_get_ksi(const union current_decrease_u * decrease, double l){
	double ans = GSL_NAN;
	switch(decrease -> id){
		case CURRENT_DECREASE_NONE:{
			ans = 1.0;
		}
		break;
		case CURRENT_DECREASE_EXPONENTIAL:{
			double eta = decrease -> exponential.eta;
			ans = exp(-eta*l);
		}
		break;
		default:
			assert(FALSE);
	}
	return ans;
}



complex current_get_spectrum(const struct current * in_current, double w, double l){
	// S_w = \ksi(l) e^{i \omega \tau(l)} \Phi_\omega
	complex S_shape_w = shape_get_spectrum( & in_current -> shape, w );
	double S_decrease_w = decrease_get_ksi( & in_current -> decrease, l);
	complex S_delay_w = cexp(I*w*delay_get_tau( & in_current -> delay, l));
	
	return S_decrease_w * S_delay_w * S_shape_w ;
}



/**
 * For given monotonically decreasing function f such that:
 * 1) f: [0,+infty) \mapsto (0,1];
 * 2) lim\limits_{x -> +\infty} f(x) = 0;
 * 3) f(0) = 1;
 * finds an argument x_{th} such that f(x_{th}) < eps.
 * Define x_1: f(x_1) = eps, x_1 - true solution.
 * The resulting x_{th} satisfies condition: x_1 < x_{th} < 1.02 * x_1.
 * It's advised that function lay within x ~ 1: e.g. if it's supposed to
 * study f(\kappa*t), kappa ~ 10^{10}, study f(x) instead of 
 * g(t) = f(\kappa*t).
 * 
 * @param f - function to study
 * @param ctx - function context
 * @param eps - treshhold criteria, eps \in (0,1)
 * @return x_{th}, discussed above
 */
static double calculate_treshhold(double (*f)(double x, const void * ctx), const void * ctx, double eps){
	double catch[2]; // a diapason such that catch[0] <= x_1 < catch[1] (see abobe)
	{
		double k = 2.0;
		double X = 1.0;
		double prevX = GSL_NAN;
		int8_t prev = 0;// {+1,-1} at last iteration f(X) {>=,<} eps; 
						// 0 - undefined, first run
		for(;;){
			int8_t this;
			double u = f(X, ctx);
			if (u >= eps){
				this = 1;
			}else{
				this = -1;
			}
			if(prev){
				if(prev != this){
					catch[0] = X;
					catch[1] = prevX;
					if(catch[0] > catch[1]){
						double tmp;
						tmp = catch[1];
						catch[1] = catch[0];
						catch[0] = tmp;
					}
					break;
				}
			} // else: this is first run
			prevX = X;
			if(this > 0){
				X = X * k;
			}else{
				X = X * (1/k);
			}
			prev = this;
		}
	}
	// We've catch[] - ed the treshhold x, catch[1] = k*catch[0].
	// We've at most (k) times overestimated it. Let's iterate a bit to make
	// the estimation closer
	//
	// ] len(catch) <= k * x_1
	// after iterations:
	// len(catch) <= k * x_1 * 2^(-nitt)  //@(k=2,nitt=7) len(catch) <= 0.016 * x_1: we won't overhelm
	int itt, nitt = 7;
	double fcatch[2];
	fcatch[0] = f(catch[0],ctx);
	fcatch[1] = f(catch[1],ctx);
	for(itt = 0; itt < nitt; itt += 1){
		double x_center, f_center;
		x_center = (catch[0]+catch[1])/2;
		f_center = f(x_center,ctx);
		uint8_t cpy;
		if( f_center >= eps ){
			 cpy = 0;
		}else{
			 cpy = 1;
		}
		catch[cpy] = x_center;
		fcatch[cpy] = f_center;
	}
	//catch diapason is 2^nitr times smaller now, let's return the upper bound as the answer
	return catch[1];
}

/**
 * Auxiliary calculator current models to be 
 * used from current_get_duration(...). 
 */
static double energy_cutoff(double x, const void * ctx){
	double ans = GSL_NAN;
	const union current_shape_u * cur = ctx;
	switch(cur -> id){
		case CURRENT_SHAPE_EXPEXP:{
			double alpha = cur -> expexp.alpha;
			double beta = cur -> expexp.beta;
			//x = alpha*t
			double t = x/alpha;
			ans = (exp(-2*alpha*t)/(2*alpha) + exp(-2*beta*t)/(2*beta) - 2*exp(-(alpha+beta)*t)/(alpha+beta)) 
				/ (1/(2*alpha) + 1/(2*beta) - 2/(alpha + beta));
		}
		break;
		case CURRENT_SHAPE_POWEXP:{
			uint8_t n = cur -> powexp.n;
			uint16_t k;
			double sum = 0;
			for(k=0;k<=2*n;k+=1){
				sum += pow((2*x), k) / gsl_sf_fact(k);
			}
			ans = sum * exp(-2*x);
		}
		break;
		default:
			assert(FALSE);
	}
	return ans;
}

double current_get_duration(const struct current * cur, double l, double eps){
	double duration_x;
	duration_x = calculate_treshhold(energy_cutoff, & cur -> shape, eps); // Signal @l=0 is tested this way for the cutoff
	
	double duration = GSL_NAN;
	switch(cur -> shape.id){
		case CURRENT_SHAPE_EXPEXP:
			duration = duration_x / cur -> shape.expexp.alpha;
		break;
		case CURRENT_SHAPE_POWEXP:
			duration = duration_x / cur -> shape.powexp.alpha;
		break;
		default:
			assert(FALSE);
	}
	
	double tau_delay;
	tau_delay = delay_get_tau(&cur -> delay, l);
	
	return duration + tau_delay; // Estimation from the above
}

const char * current_validate(const struct current * cur){
	////
	// Checking shape
	////
	const union current_shape_u * shape = & cur -> shape;
	switch(shape -> id){
		case CURRENT_SHAPE_EXPEXP:{
			double alpha = shape -> expexp.alpha, beta = shape -> expexp.beta;
			if(!(alpha < beta) || !(alpha>0) || !(beta>0)){
				return "alpha must be less than beta, and both positive";
			}
		}
		break;
		case CURRENT_SHAPE_POWEXP:{
			double alpha = shape -> powexp.alpha;
			uint8_t n = shape -> powexp.n;
			if ( alpha <= 0 ){
				return "alpha must be positive";
			}
			if( n == 0 || n > 10 ){
				return "n must be positive, reasonably small integer";
			}
		}
		break;
		default:
			return "Unknown shape type";
	}
	////
	// Checking delay
	////
	const union current_delay_u * delay = & cur -> delay;
	switch(delay -> id){
		case CURRENT_DELAY_LINEAR:{
			double v0 = delay -> linear.v0;
			if( v0 <= 0 ){
				return "v0 must be positive";
			}
		}
		break;
		case CURRENT_DELAY_EXPONENTIAL:{
			double v0 = delay -> exponential.v0;
			double gamma = delay -> exponential.gamma;
			if(v0 <= 0){
				return "v0 must be positive";
			}
			if(gamma <= 0){
				return "gamma must be positive";
			}
		}
		break;
		default:
			return "Unknown delay type";
	}
	////
	// Checking decrease
	////
	const union current_decrease_u * decrease = & cur -> decrease;
	switch(decrease -> id){
		case CURRENT_DECREASE_NONE:{
			// Nothing to check :)
		}
		break;
		case CURRENT_DECREASE_EXPONENTIAL:{
			double eta = decrease -> exponential.eta;
			if(eta < 0){
				return "Eta must be non-negative";
			}
		}
		break;
		default:
			assert(FALSE);
	}
	
	// OK, returning NULL
	const char * ans = NULL;
	return ans;
}




/**
 * Calculates the values to find treshhold of for some of the currents.
 * @param w - 2*pi*frequency
 * @param ctx - current object pointer
 * @return value to find treshhold of
 */
static double get_max_w_func(double w, const void * ctx){
	const union current_shape_u * sha = ctx;
	double ans = GSL_NAN;
	switch(sha -> id){
		case CURRENT_SHAPE_EXPEXP:{
			double alpha = sha -> expexp.alpha;
			double beta = sha -> expexp.beta;
			ans = alpha*beta / sqrt(  (pow(beta,2) + pow(w,2)) * (pow(alpha,2) + pow(w,2))  );
		}
		break;
		case CURRENT_SHAPE_POWEXP:{
			double alpha = sha -> powexp.alpha;
			uint8_t n = sha -> powexp.n;
			ans = pow(alpha, n+1)/pow(pow(alpha,2) + pow(w,2),(n+1)/2.0);
		}
		break;
		default:
			assert(FALSE);
	}
	return ans;
}


double current_get_max_w(const struct current * cur, double eps){
	double ans;
	ans = calculate_treshhold(get_max_w_func, & cur -> shape, eps);
	return ans;
}


static void current_shape_mirror(const union current_shape_u * src, union current_shape_u * dst){
	* dst = * src ; // Binary copy
	switch(src -> id){
		case CURRENT_SHAPE_EXPEXP:{
			dst -> expexp.Q = (-1) * src -> expexp.Q;
		}
		break;
		case CURRENT_SHAPE_POWEXP:{
			dst -> powexp.Q = (-1) * src -> powexp.Q;
		}
		break;
		default:
			assert(FALSE);
	}
}


void current_mirror(const struct current * src, struct current * refl){
	refl -> delay = src -> delay;//Copy binary data!
	refl -> decrease = src -> decrease;
	current_shape_mirror(& src -> shape, & refl -> shape);
}





